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Many-body ground states can be prepared via unitary evolution in cold atomic systems. Given the initial state 
and a fixed time for the evolution, how close can we get to a desired ground state if we can tune the Hamiltonian 
in time? Here we study this optimal control problem focusing on Luttinger liquids with tunable interactions. 
We show that the optimal protocol can be obtained by simulated annealing. We find that the optimal interaction 
strength of the Luttinger liquid can have a nonmonotonic time dependence. Moreover, the system exhibits a 
marked transition when the ratio t/L of the preparation time to the system size exceeds a critical value. In this 
regime, the optimal protocols can prepare the states with almost perfect accuracy. The optimal protocols are 
robust against dynamical noise. 



Loading cold atomic gases into optical lattices provides an 
opportunity to study the nonequilibrium properties of quan- 
tum matter in thermally isolated and highly tunable environ- 
ments [1,2]. The central object of such studies is a many-body 
quantum state which undergoes unitary evolution generated 
by a time-dependent local Hamiltonian. 

In an important class of problems, we are interested in us- 
ing the unitary evolution to transform an initial state that is 
the ground state of a local Hamiltonian to the ground state 
of another Hamiltonian. Such problems are relevant for the 
preparation of states in regimes where direct equilibration is 
difficult [3-5]. Ground states preparation is the key to simu- 
lating many-body model Hamiltonians, such as the Hubbard 
model, which can be realized with cold atoms [6, 7]. 

If one had infinite time to wait, according to the adiabatic 
theorem of quantum mechanics, the unitary transformation 
can be done with arbitrary accuracy in any finite system. Ex- 
trinsic losses and quantum decoherence, however, set an upper 
bound on the practical time to carry out the process. In any fi- 
nite time t, nonadiabatic effects are unavoidable [8]. These 
effects are most severe in the absence of an energy gap [9]. In 
this paper we focus precisely on dynamics within the gapless 
phase by studying Luttinger liquids. 

Let us start by casting the question of finding the opti- 
mal dynamical protocol in a generic way. Assume we have 
a local Hamiltonian H({g}) = YldLi 9i ®i where the OiS 
are local operators and the <^s are coupling constants that, 
within a given range, can be tuned to any value as a func- 
tion of time. We would like to transform which is 
the ground state of H({gi}) to |v&2), the ground state of 
H({g2}) in a given time r. How close can the final state 
|*(r)) = Texp[-i /J" (ft' J ff({ 3 (t')})]|1 , i) be to the desired 
ground state | SI/2) ? Here Trepresents time-ordering. 

The meaning of closeness above depends on the measure 
used. There are several popular measures like the excess en- 
ergy for example [10]. Here we use the wave function overlap 



F{{g(t)}] = \(*(T)\*2)\ 2 



(i) 



The problem is then reduced to finding the time-dependent 
{g{t)} that maximizes the functional above. This interesting 
question in quantum dynamics [8] is in fact a typical problem 
in optimal control theory as noted in Ref. [11, 12]. Let us 



emphasize that we are concerned only with the final state and 
maintaining adiabaticity during the evolution, as for example 
in Ref. [13], is not a constraint. Moreover, we are restricted to 
local Hamiltonians allowed by the experimental setup. 

The focus of this work is the Luttinger model for interacting 
fermions in one space dimension, whose Hamiltonian can be 
written in momentum space as follows 



q>0 ^ 



(2) 



where $ 9 are bosonic fields and U q their conjugate momenta. 
The parameters u and K are respectively the velocity of the 
charge carriers and the Luttinger parameter. We consider a 
fixed number of particles, focusing on the half-filled case. 
The zero mode, which is responsible for changing the parti- 
cle number sector in the bosonized description, is therefore 
excluded in the above expression. Assuming we have an odd 
number of sites L in the system, the momenta q are given by 
9 = 27rf forn = l--- , 

This Luttinger Hamiltonian is the low-energy effective the- 
ory for models of interacting fermions on a one-dimensional 
lattice, in particular the ID Hubbard model 



-c]c j+1 + h.c. + V{nj - - 



(3) 



Since spin degrees of freedom are not essential for our discus- 
sion, here we focus on spinless fermions. With the hopping 
amplitude set to unity as in Eq. (3), the Luttinger parameters 
u and K are related to V via the Bethe ansatz (see Ref. [14] 
for example). 

We will consider trajectories for u and K that are 
parametrized by a time-dependent V(t): u — u{V) and 
K = K{V). The strength of Hubbard-type interactions can 
be tuned in optical lattices both by manipulating the opti- 
cal potential with lasers [15-17] and through Feshbach res- 
onances controlled with magnetic fields [17, 18]. 

Notice that while the relation between the Luttinger liquid 
(LL) and the ID Hubbard model holds at low energies, the 
results we find for the optimum dynamical protocol in the for- 
mer should be applicable to the latter when the total momen- 
tum qn q in each harmonic oscillator mode in Eq. (2), where 



2 



n q is the occupation number, is small compared with ir/2 [19]. 
We shall check this condition a posteriori. 

At half filling, the gapless LL description holds for -2 < 
V < 2. At V = 2, a charge density wave (CDW) gap 
opens up. An optimal power-law protocol for bringing the 
system from the gapped phase to the critical point was found 
in Ref. [20] using adiabatic perturbation theory [21]. Here we 
consider the problem of transforming the system initially at 
the CDW phase transition critical point to a point deep within 
the gapless LL phase. 

Let us now assume we are initially in the ground state at 
the critical point (Vi = 2, K\ = \). We would like to find 
the time-dependent interaction strength —2 < V(t) < 2 for 
< t < t that yields the maximum overlap with the ground 
state of the Hamiltonian Eq. (3) with V = Vi at time r. 

We proceed by expressing H q = —i d<s, q in Eq. (2). The 
time-dependent many-body wave function of the system for a 
protocol V{t) (and consequently u{t) and K{t)) can then be 
written in the |{$ q }) basis as 

*({*,» - n [^otf ^ *) <m)] (4) 

q>0 

where Sft and 9 indicate the real and imaginary part and 
^ q ((j), t) is the solution of the following Schrodinger equation 



id t - u(t) 



1 



K(t) 



= 



with appropriate initial conditions. 

Up to an unimportant overall phase, the solution of the 
above differential equation is 



2q 



[5ftz 9 (t)]iexp(- g z 9 (i)0 2 ) (5) 



where z q (t) is a complex-valued function that satisfies the fol- 
lowing Riccati equation 



iz q {t) = q^ ) [zl{t)-a\t)] 



(6) 



with a(t) = 1/K(t) and the initial condition z q (0) = 

To perform the optimization for the many-body system, we 
discretize time and approximate a general V(t) by a piece- 
wise constant function over the interval [0 . . . r]. This allows 
us to write the final overlap, which is a functional of V(t), as 
a multi-variable function that can be maximized numerically. 
An unbiased optimal protocol can be found by increasing the 
number of discretization points. 

Note that our approach to computing the overlap could also 
be useful for a variety of other calculations (see [22, 23] for 
example) in the quench dynamics of Luttinger liquids. Let 
us assume a sequence Vj with j = 1 . . . N such that with 
At = t/N, V(t) = Vj for (j - l)At < t < jAt. We then 
get two corresponding sequences Uj and aj. If z q = z q (jAt), 



we obtain the following recursion relation for z q by solving 
Eq. (6) for time-independent u and a, 



i otj tan 



q Uj At + arctan —i- 



7? : -i 



(7) 



Our focus here is finding an optimal protocol but Eq. (7) 
above is of interest in its own right since it gives an exact solu- 
tion of the nonequilibrium wave function for any sequence of 
sudden quenches in the Luttinger model. Notice that knowing 
z q (t) for all modes determines the many-body wave function. 

Recursively solving the above relation Eq. (7) yields 



z q (r) 



for any given piece-wise constant interaction 



strength. Defining a 2 
be written as 



J-(Vi,...V 2 v)=exp 



1/K 2 , the overlap Eq. (1) can then 



q>0 



4 a 2 



ftja(r) 
"2 + z q (r) 



(8) 

The overlap above is written as a multi-variable function 
of {Vj}. To find the optimal protocol, we minimize the cost 
function £({Vj}) = — In ^({Vj}) with respect to the config- 
uration {Vj}. This can be done by Monte-Carlo (MC) meth- 
ods. We perform simulated annealing calculations with ki- 
netic moves consisting of random small changes in randomly 
chosen VjS. 

We compare the results of the optimal protocol against two 
additional calculations. We consider the one-parameter varia- 
tional protocol V(t) = Vi + (V 2 - Vi) (t/r) r and calculate the 
final overlap for the linear (r = 1) as well as for the the best 
power-law protocol (r = r m ; n with d r £{r) = 0). The opti- 
mal protocol found by MC simulations performs significantly 
better than both of the above. 

In Fig. 1, we show the final £ for the optimal protocol ob- 
tained by an unbiased MC simulation as well as for the linear 
and the best power-law protocols for V 2 = —1-5 and several 
system sizes. When t/L becomes larger than a critical value, 
the final £ obtained by MC optimization exhibits a qualita- 
tive change of behavior and shoots down by several orders of 
magnitude. 

The value of the cost function obtained by MC simulations 
is generically only an upper bound on the actual minimum. 
For r < t c oc L, different annealing histories lead to a unique 
protocol suggesting that the the minimum found is likely the 
global minimum. For r > t c however, we never converge 
to a unique protocol indicating that we have only found a lo- 
cal minimum. In this regime, finding the global minimum 
becomes exceedingly difficult, particularly for larger systems 
due to the increased computational complexity. Nevertheless, 
notice that the minima we find have £ (r) very close to zero 
and can prepare the new ground state almost exactly. 

Although it is not clear from the finite-size scaling of the 
obtained £{t) that there is a transition in the thermodynamic 
limit, the drastic change of the system behavior strongly sug- 
gests the presence of a true transition for the global minimum. 
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FIG. 1. The final £ obtained by Monte-Carlo (MC) simulations 
compared with the linear and best power-law protocols. The final 
£ plunges by several orders of magnitude for t/L > t c /L w | 
(marked with a star). Note that in this regime, the data merely repre- 
sents an upper bound on £(t) and the actual cost function could be 
even smaller. 
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FIG. 2. The interaction strength V(t) as a function of time t for 
the optimal protocol obtained by MC simulations for L = 129, 
t = 8 and V% = —1.5. Here, the optimal protocol converges to a 
smooth oscillatory function as we increase the number of discretiza- 
tion points N. 

Note that the local minima we find are expected to be closer 
to the actual global minimum for smaller systems. 

An interesting feature of the optimal protocols is that, as 
seen in Fig. 2, the interaction strength can be an oscillatory 
function of time. Note that increasing the number of dis- 
cretization points leads to convergence to smooth but non- 
monotonic protocols in the r < r c regime. For large enough 
systems, the period of the oscillations does not have a strong 
dependence on the system size L or the preparation time r. 
With no other time scale left, we then conjecture that the os- 
cillations must be a consequence of the short-distance length 
scale, i.e. the lattice spacing which is set to unity in our prob- 
lem. As seen in Fig. 3, the period of the oscillations decreases 
as Vi becomes larger. This observation is consistent with a 
short-distance cut-off controlling the oscillations as the veloc- 
ity u(V) is an increasing function of the interaction strength 
V. 

We now speculate a possible scenario for the nature of the 
transition above. Let us consider a single mode q. We ex- 
pect to be able to find an exact solution to z q (r) = a 2 by 
using a two-parameter variational protocol and solving two 
equations (real and imaginary part of z q ) and two unknowns. 
Genetically however, this system of equations does not ad- 
mit a solution. For r = 0, for example, z q cannot change 
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FIG. 3. The optimal protocol interaction strength V(t) for two dif- 
ferent values of Vu = —1.0 and V2 = —0.5. 

from z q (0) = a.\. By rescaling time in Eq. (6), we notice 
that the existence of a solution depends only on the quantity 
rq. We have checked numerically that with a two-parameter 
piece-wise constant protocol, we are able to obtain solutions 
to z q (r) = a 2 with \a 2 — ol\\ ~ 0(1) only when qr > 0(1). 

We conjecture that even with an infinite number of varia- 
tional parameters, a minimum time of 0(l/q) is still required. 
Circumstantial support for this conjecture comes from the MC 
simulations with piece-wise constant protocols where, as seen 
in Fig. 2 for t < r c , there is convergence in the shape of the 
optimal protocols as we increase N. Notice that although we 
have more variation parameters Vj for larger N, each of them 
acts for a shorter time. When qr is large enough such that, 
with two variational parameters, z q (r) = a 2 has a solution, 
we have an infinite number of solutions for N > 2. 

For the many-mode problem, we wish to have z q (r) = a 2 
for all the modes. From the discussion of the single-mode case 
above, we find that, for times larger than some r oc L, an infi- 
nite number of solutions exists for the slowest mode (q = 
and consequently for every individual mode q. Since the sys- 
tem of equations z q (r) = a,%, Vg is under-determined when 
the number of variational parameters is larger than L — 1, it 
seems plausible that for r > r c oc L, one can simultaneously 
satisfy z q {r) = a 2 for all modes. We have not, however, been 
able to explicitly find such exact solutions for the many-mode 
problem and thus cannot rule out the possibility that the tran- 
sition at t > t c (x L happens because \z q (r) — a 2 \ for all 
modes can become infinitesimally small rather than exactly 
zero. 

Before proceeding, let us check the applicability of the re- 
sults obtained for the Luttinger model Eq. (2) to the Hubbard 
model Eq. (3). We need to calculate the total mode momentum 
qn q . Writing the occupation number as n q = (e q ) /2e q — 1/2 
where e° is the ground-state energy of a mode, we obtain 

= (*(*) MOI 2 + jfa) ~ §■ We find that 

when the final overlap is large, the evolution does not typi- 
cally excite electrons too far away from the Fermi surface. For 
example, the protocol for L = 129 and t = 16 in Fig. 2 has 
max g jt [q n q(t)] — 0.302 which is indeed in the linear regime. 

From the experimental point of view, the extremely small 
£ (t) obtained in the regime r > t c can be very useful for the 
preparation of many-body states. In practice however, there 
are always inaccuracies in the experimental implementation 
of a prescribed protocol. We check the robustness of these 
protocols against random perturbations 6V(t) taken from a 
uniform distribution [— W/2 . . . W/2] for each segment of the 
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FIG. 4. Left: the noise-averaged £(t) as a function of W 
for a random noise 5V(t) taken from a uniform distribution 
[-W/2 . . . W/2\. Right: final cost £ (r) as a function of 8V\ where 
the initial wave function is the ground state for interaction strength 
Vi — SVi. The plots are for L = 65, r = 16 and Vi = —1.5. 



piece-wise constant V(t), We also check the effect of errors 
in the initial wave function by applying the protocol to the 
ground state for V\ — SVi instead of V\. 

In Fig. 4, we show £(t) as a function of SVi in addition to 
the noise-averaged £ (r) as a function of W for the MC opti- 
mized, best power-law and linear protocols and for L = 65, 
r = 16 and V% — —1.5. As seen in the figure, even for large 
noise of the order a few percent of the bandwidth, the MC op- 
timized protocols yield much smaller £ (t) than the power-law 
and linear V(t). These studies indicate that, practically, using 
the optimal protocols is advantageous even in the presence of 
experimental errors and inaccuracies. 

To better understand the effects of this dynamical noise on 
the behavior of thermally isolated Luttinger liquids, we con- 
sider a simple case. We assume we are initially in the ground 
state for interaction strength Vq (with corresponding «o and 
uq) and V(t) — Vq + SV(t) for a time r where 5V(t) is a 
random noise of strength W. For W — 0, the wave function 
remains in the ground state and £ (r) = 0. For small W, we 
can linearize Eq. (6) and write iSz q = 2quQ (8z q — 8a) where 
8z q = z q — a and 8a = a — a Q . We then obtain 

8z q (T) = 2iquo [ dte 2tqUo{t ~ T) 8a{t). (9) 
Jo 

Using Eq. (1), we get £{r) = ^ E 9 >ol^WI 2 + 
0(8z 3 ) which assuming uncorrected noise (8a(t)8a(t')) oc 
W 2 8(t — t 1 ) and using Eq. (9) above leads to £(t) oc tW 2 . 
Note that since r c oc L, strong noise limits the size of the 
systems which can be accurately prepared. Also notice that 
|<5z,j(t)| 2 oc q 2 and in contrast to coupling to a thermal bath, 
the dynamical noise discussed above creates more excitations 
in modes with higher momentum. The effect of noise on a 
time-dependent optimal protocol seen in Fig. 4 is mathemat- 
ically more complicated but essentially similar to this simple 
case. The noise-averaged £(t) grows as a power law with W 
with an exponent that depends on V\ and V%. 

In summary, we used simulated annealing to address an 
outstanding problem in the nonequilibrium dynamics of in- 
teracting quantum systems, namely finding unbiased optimal 
protocols for unitary preparation of strongly correlated states. 
We focused on transforming states in the LL phase of inter- 
acting fermions with tunable interaction strength. Unbiased 



optimization over all ramp shapes yields nonmonotonic opti- 
mal protocols. The behavior of the system exhibits a marked 
transition for r > r c oc L where the states can be transformed 
with almost perfect accuracy. The optimal protocols are robust 
against noise, which makes them of practical experimental im- 
portance for the preparation of states in cold atomic systems. 
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problem is studied for small ramps in infinite dimensions with 
dynamical mean field theory. This work was supported by the 
DOE Grant DE-FG02-06ER463 16. 
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Abstract 

The trajectories of z q (t) on the complex z plane are discussed for the 
evolution with optimal protocols. 

Some intuition may be gained regarding the appearance of oscillations in 
optimal protocols discussed in the main text by considering the trajectories of 
z q (t) on the complex z plane. The evolution of the system is governed by a set 
of nonlinear differential equations (Eq. (6) of the main text) for many modes. 
We would like controls a(t) and u(t) such that for all modes q, the complex 
function z q (t) moves from the initial value z q (Q) = otx to a desired final value 
ai in a time r. Now the quantity controlling this evolution for each mode is qr 
as discussed in the main text. It is evident from Eq. (6) of the main text that 
z q moves faster for modes with higher momenta. To get to the same point ct2, 
we find that z q actually winds around in the complex plane as seen in Fig. 1. 
The appearance of these windings, i.e oscillations in the imaginary part of z q , 
is rather intuitive. One may guess that an oscillatory control is required to 
generate the winding of the complex z q (t) . 
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Figure 1: The values of z q (t) for three different modes during the evolution 
with the MC optimized protocol with L = 65, r — 16 and Vi ~ —1.5. As the 
momentum q increases the modes wind around in the complex plane to reach 
the point a-i- 



